Tidal interaction in compact binaries: a post-Newtonian afRne framework 
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We develop a semi-analytical approach, based on the post-Newtonian expansion and on the afBne 
approximation, to model the tidal deformation of neutron stars in the coalescence of black hole- 
neutron star or neutron star-neutron star binaries. Our equations describe, in a unified framework, 
both the system orbital evolution, and the neutron star deformations. These are driven by the 
tidal tensor, which we expand at 1/c^ post-Newtonian order, including spin terms. We test the 
theoretical framework by simulating black hole-neutron star coalescence up to the onset of mass 
shedding, which we determine by comparing the shape of the star with the Roche lobe. We validate 
our approach by comparing our results with those of fully relativistic, numerical simulations. 

PACS numbers: 04.25.Nx, 04.30.Dg, 04.25.dk 



I. INTRODUCTION 

Coalescing binaries composed of neutron stars (NS) 
and/or black holes (BH), are among the most promising 
sources of gravitational waves (GWs) to be detected by 

Savitational wave interferometers like Virgo and LIGO 
. These systems are also interesting since they are 
thought to be related to short gamma-ray burst . 

The process of coalescence has been studied mainly 
using post-Newtonian (PN) and fully relativistic, numer- 
ical simulations. PN expansion Q has the advantage 
of providing a semi-analytic description of the evolution 
of the system, but it is poorly convergent in the strong 
field limit; therefore it is appropriate to study the inspiral 
phase only. These limitations have been overcome by PN 
resummed formulations like BOB but these formula- 
tions are not able to describe the dynamical features of 
the stellar deformation. Moreover, in the standard PN 
expansion the compact objects are treated as pointlike up 
to the 4.5 (included) post-Newtonian order. Finite size 
effects are, formally, of order 5PN; however their contri- 
bution is larger than what a naive counting of PN orders 
may suggest Q. Tidal deformations have recently been 
included in the PN framework through the "Love num- 
ber" approach 01 , which assumes that the tidal field is 
proportional to the quadrupole momentum (see below). 
The effects of the tidal deformation on the orbital motion 
have been studied in [1, @ . 

Fully relativistic codes are the most powerful tool to 
investigate the latest phases of the inspiral and merger 
(see [lOl for a review on the subject). They are, how- 
ever, not exempt from drawbacks: their computational 
cost is high, therefore the parameter space cannot be ex- 
plored at large; furthermore, initial data solvers are still 
unable to provide accurate initial data for binaries with 
non-aligned spins, and may introduce spurious numerical 
effects which, if not appropriately cured, affect the sub- 
sequent evolution of the system. These problems are of 
particular relevance in BH-NS binaries, where the lack of 
symmetry makes more difficult to follow the entire pro- 
cess of coalescence by fully relativistic simulations. For 
these reasons, the process has been studied in the liter- 



ature using some simplifying assumptions, or for a re- 
stricted set of parameters. For instance, in (lT| - [l3 | the 
inspiral is modeled as a sequence of quasi-equilibrium 
circular orbits with decreasing radius; in [TTI.[T2j the pro- 
cess is studied by fully relativistic simulations, whereas 
[H, [i3l use the affine approach (see below) . In [l5|, [l^ 
Einstein's equations are evolved assuming that the black 
hole is non rotating, and for large values of the mass ra- 
tio q ~ Mbh/Mns, whereas in p7l - l20j q takes values 
g < 5; in [2l| the black hole is assumed to rotate 
with spin parallel to the orbital angular momentum, and 
different values have been considered. For a recent re- 
view on fully relativistic simulations of BH-NS binaries 
see [13 (the literature on NS-NS coalescing binaries is 
much more extended, and we do not report it here). 

In this paper we develop a semi-analytic approach to 
study BH-NS and NS-NS coalescence, by merging two 
different frameworks: the PN approach, which accurately 
describes the system orbital motion, and the affine model 
[Tsl . [1^ [23l - l26j , which describes the stellar deformations 
induced by the tidal field. To this aim, we compute the 
tidal tensor associated to the PN metric of a two-body 
system, defined in terms of the PN Riemann tensor and 
of the local tetrad of the deformed body 

C(,)(j) = Ra(s^sefo)e1,fJoffj) , (1) 

up to 0(l/c^). This tensor was derived with a different 
approach in [27h29| up to 0(l/c^); our expression coin- 
cides with that of [23 and also includes 0{l/c^) terms, 
associated to the spins of the compact objects. 

In the affine model the NS is described as a deformable 
ellipsoid, subject to its self-gravity, to internal pressure 
forces and to the tidal field of the companion. In the orig- 
inal formulation of this approach, the NS structure was 
considered at a Newtonian order, assuming a polytropic 
equation of state (EO S) J23l - [26j . A first improvement 
was introduced in [13l . uM , where general relativity was 
taken into account in the description of the stellar struc- 
ture, and non-polytropic EOSs were considered. This 
approach was used to study quasi-equilibrium configura- 
tion sequences of BH-NS systems, in order to estimate 
the critical distance at which the NS is disrupted by 
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the tidal interaction [ij], to determine the correspond- 
ing cut-off frequency in the emitted gravitational wave 
signal [3, [131, and to estimate the mass of the torus 



which forms after the NS is disrupted [31 1. 

In this paper we further improve the afhne model intro- 
ducing a more accurate description of the orbital motion 
and of the tidal interaction. Our approach differs from 
existing work on NS tidal deformation in compact bina- 
ries in the following aspects. 



tional to the tidal tensor: 



(2) 



• In [13|,llJ,|25|,|26|,|3lll32l the affine model was used 
assuming that the NS follows a timelike geodesic 
of Kerr's spacetime; this approximation fails when 
the mass ratio q is low. In addition, time-dilation 
factors were neglected. 

Furthermore, most works employing the afhne 
model [H, [ij, HE do not evolve the dynami- 
cal equations of stellar deformation. Rather, they 
find, at each value of the orbital radius, the cor- 
responding stationary configuration describing the 
deformed star. In [U 113 the dynamical equations 
were solved; however, while the orbital evolution 
was described in PN coordinates, the BH tidal field 
was expressed in the Boyer-Lindquist coordinates 
of the Kerr metric describing a single BH; neglect- 
ing the difference between these coordinate systems 
yields a loss of accuracy of the model. 

These problems are solved in the fully consistent 
approach presented in this paper, where the BH-NS 
or NS-NS systems are described by a two-body PN 
metric, which holds for any value of the mass-ratio 
q. The tidal tensor itself is expressed in the PN 
coordinates, and the proper time of each compact 
body is expressed in terms of the PN time coordi- 
nate, through the appropriate Lorentz factor. Our 
approach is valid up to the onset of mass shedding, 
which occurs when the deformed star crosses the 
Roche lobe; after that, it can no longer be applied, 
since the assumption that the star is a deformed 
ellipsoid is significantly violated. We describe the 
orbital motion of the compact objects using the 3.5 
PN equations for pointlike objects, with next-to- 
leading-order tidal corrections; the NS tidal defor- 
mation is driven by the tidal tensor of the 3 PN 
metric. The dynamical equations are a system of 
(non-linear) ordinary differential equations in time. 

• NS tidal deformations have also been studied in a 
series of paper [1, 0, [H-lHl , where the deformation 
properties have been encoded in a set of numbers, 
the Love numbers, which relate the quadrupole ten- 
sor (or, more generally, the multipole moments of 
the star) to the tidal tensor. This approach is 
grounded on the adiabatic approximation, i.e., on 
the assumption that the orbital evolution timescale 
is much larger than the timescale needed for the 
star to set into a stationary configuration. In this 
approximation, the quadrupole tensor is propor- 



with A constant. The Love number A can be com- 
puted by studying the response of a single star to 
an external tidal tensor 0, 0, 113, HI] ■ This model 
has been employed to determine the effect of tidal 
deformation on the orbital motion of a NS in a bi- 
nary system d, [H, [H [H] . 

We also compute the Love number A (see Sec- 
tion imp , without assuming the adiabatic approxi- 
mation: the stellar deformation is found by solving 
dynamical equations. 

To test the accuracy of our approach, we compare the 
results with the existing literature on BH-NS binaries. 
As a preliminary check, we verify that our PN descrip- 
tion of the orbital motion accurately reproduces the fully 
relativistic results [H, [10, [13 • Then, we verify that the 
onset of mass shedding we determine, is consistent with 
the results of fully relativistic simulations [l^, [l^, [STj . 
Finally, we check that the stellar deformations predicted 
by our model are consistent with existing computations 
of the Love number ^ [3, [H, [Hi- 

This paper focuses mainly on the theoretical frame- 
work, and on its validation by comparison with the ex- 
isting literature, where available. The tool we develop 
will be used in future works to study the dynamics of 
compact binaries. 

The plan of the paper is the following. In Sec HI] we de- 
scribe the model. In Sec. Illll we assess the validity of our 
approach, by comparing the results with the literature. 
In Sec. IIVI we draw the conclusions. 



II. THE MODEL 

We use notations and conventions introduced in [l3| . 
where the affine model (partially improved with respect 
to the original formulation [l^, [3] ) is widely discussed. 
In the following mi,m2 are the masses of the two com- 
pact objects; we shall consider the tidal deformation of 
the NS with mass mi and radius Rns] the companion, 
with mass m2, can either be a BH or a NS. Furthermore, 
we define m = mi -|- m2 and u = ^/m ^ mim2/m^, and 
the mass ratio q = mijmx. 



A. Improved afflne model 

The basic assumption of the affine model is that the 
NS, deformed by the tidal field, maintain an ellipsoidal 
shape; it is an 5'- type Riemann ellipsoid, i.e., its spin 
and vorticity are parallel, and their ratio is constant [38| . 
The deformation equations are written in the star princi- 
pal frame, i.e., the frame comoving with the star, whose 
axes coincide with the ellipsoid principal axes. In what 
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follows, ai is the axis which points toward the compan- 
ion; 02 and 03 are the axes orthogonal to ai, with 02 
lying in the orbital plane; the indices 1,2,3 label the cor- 
responding directions. Surfaces of constant density inside 
the star form self-similar ellipsoids and the velocity of a 
fluid element is a linear function of the coordinates 
in the principal frame. Under these assumptions, the in- 
finite degrees of freedom of the stellar fluid motion can 
be reduced to five [231 - 125} , and are associated to dynam- 
ical variables governed by a set of non-linear differential 
equations, which describe the evolution of the stellar de- 
formation. These variables are the three principal axes of 
the ellipsoid ai {i = 1, 2, 3) and two angles, tp, A, defined 
as 



dr 



dX . 



(3) 



where r is the NS proper time, and VL is the ellipsoid an- 
gular velocity, measured in the parallel transported frame 
associated with the star center of mass, i/' is the angle 
between the principal frame and the parallel transported 
frame. A is defined as follows: 



A 



0102 



as 



■c 



(4) 



where C, is the vorticity along the axis in the principal 
frame. The NS internal dynamics is described in terms 
of the Lagrangian 



Ci=Ti~U-V , 



(5) 



where T/ is the fluid kinetic energy, U is the internal 
energy and V is the star self-gravity. In the original ap- 
proach introduced by Carter and Luminet, these are de- 
fined in a Newtonian framework 



Ti = \ I v'dMs , 



U 



V = - 



-dMp 



G f dMsdM'g 



= / dMerdr^ 



Newt 7 



(6) 
(7) 
(8) 



where Alg is the baryonic mass, p the mass density, e the 
Newtonian energy density, ^Newt the gravitational po- 
tential; all these quantities are solutions of the Newtonian 
equations of stellar structure. Furthermore, dMs = pd^x 
and V satisfies the virial theorem, which states that, in 
the spherical configuration, V = —311, where 



n = / -dMe = 
P 



pd^x . 



(9) 



The variation of the Lagrangian ([S]) gives the equations 
of motion for the five dynamical variables a^, ip, X. 

In [2^, it was shown that, under the affine hypoth- 
esis, the integrals in Eqns. (jBHEl) and their variations, can 
be expressed in terms of integrals on the fiuid variables 
computed for the spherical configuration (a^ = R^s), 



and of functions of the dynamical variables. With this 
simplification, the equations of motion can easily be 
found. In the following, a superscript hat will denote 
quantities computed for the spherical star. T and V ex- 
pressed in terms of the "hatted" quantities and of the 
dynamical variables are 



Tr = 



1 



^ 2 
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1 M 
2 

1 



daj 



^A 

a2 



M 
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0-2 



il A 

ai 

da 



where A4 is the scalar quadrupole moment 
1 



M 



sph 



(10) 
(11) 

(12) 



(the subscript sph means that the integration is per- 
formed on the spherical star) and V is the self-gravity 
potential of the spherical star. 

The procedure to make explicit the dependence of the 
internal energy U on the dynamical variables is more sub- 
tle. The internal energy variation dU can be written as 



n 



dU = — dai . 



(13) 



The pressure integral 11 given by Eq. (jH]) can not be fac- 
torized in a spherical integral and a function of the axes; 
however, it can be expressed as 



n 



aia2a3 



p{p)d^x 



(14) 



sph 



where p{p) is the fluid equation of state, and p is the 
rescaled mass density 
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(15) 



with p mass density in the spherical configuration. When 
Oi = i?Ar5, the pressure integral 11 reduces to the spheri- 
cal pressure integral 11. 

A first improvement to this approach was introduced in 
[l3| , where the Newtonian description of the NS equilib- 
rium configuration was replaced by the relativistic equa- 
tions of stellar structure (TOV) 



dr 



dr,,P = —G 



{e + p/c'^){ms + ATiprl/c^) 
rs{rs - 2Gms/c?) 



(16) 



Here e is the relativistic mass-energy density in the 
spherical configuration, is the radial coordinate in a 
Schwarzschild frame associated to the non rotating NS, 
and nia (r,, ) is the gravitational mass enclosed in a sphere 
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of radius rg. We remark that this is a major change, 
since the relativistic radius is smaller than the Newto- 
nian radius by ~ 10% — 20%. Wc also remark that the 
Schwarzschild coordinate is different from the radial 
coordinate in the Newtonian frame r = The 
self-gravity integral V was changed accordingly, as 



V= / dAlBrsdr^^TOV 
' sph 



(17) 



where $Toy is an effective relativistic gravitational po- 
tential of the spherical star, defined in terms of the TOV 
equations as follows 



TOV 



^ (e + p/c^)(m, + 47rr3p/c^) ^^^^ 



rs{rs - 2Gms/c'^) 
With this definition the virial theorem 



Tsdr^^TOvdMB 



sph 



-3 / ^dMs 

'sph P 



(19) 



is still satisfied, with p solution of the TOV equations 
(jl6p and p baryon mass density. As shown in Section 
III El some terms in the dynamical equations cancel in 
the spherical limit, only if the virial theorem is satisfied. 
A non-exact cancellation of these terms would lead to 
strong instabilities. 

A further improvement, which we introduce in this pa- 
per, consists in a careful treatment of the coordinate 
frames. To describe the integrals in the spherical con- 
figuration, the relevant coordinate systems are: (i) the 
Schwarzschild frame, with radial coordinate r^, in which 
the TOV equations are expressed; (ii) the Newtonian 
frame for a spherical star {a;*}, which we now replace with 
the corresponding 1 PN post-Newtonian frame with 



isotropic radial coordinate r 
a single star) 



^-(x^) and metric (for 



1 



2V 



2V^ 



dt^ 



2V 

72" 



Sijdx^dx-' 



(20) 



where V{r) = G /~ ^^^y^dr'. Following [H the trans- 
formation between the post-Newtonian isotropic radial 
coordinate and the Schwarzschild coordinate inside the 
star is given by 



1 - 



(21) 



The scalar quantity dAfg can be expressed, in the 
Schwarzschild frame, in terms of the corresponding spa- 
tial three- metric 'jHf^^ : 



dMs = Ryjlschwd X 

2 / 1 , Gms{rs) 



prs 



1 



dr, sin 



(22) 



The integrand in the quadrupole moment (|12p is ex- 
pressed in the post-Newtonian coordinates, i.e., it is 
= r1{l - 2V{rs)/c^). The integrals V, n, M then take 
the form 



V = -3n 

n = 47r 
M 



An /"^"^ , 



2V{r,) , Gm,(r,) ^ 4 



r^drs 



(23) 



B. The post-Newtonian metric 

To derive the equations describing the orbital mo- 
tion of the binary and the tidal tensor, we shall use a 
3 PN metric written in harmonic coordinates ({a;^ = 
ct,x,y,z}): 





= -1 
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X + v,v, + - 
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32 

78 



VX 



16 
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f 0(9) 
4 



vv,y,. 



48 



0(10) 



16 



,6 \ ^y' 



0(8) 



(24) 

(25) 
(26) 



where the potentials V,Vi, X ,Wij , Ri,Yi,Zij , are defined 
in terms of retarded integrals over the source densities 
[40I l4l| . We stress that the potential V appearing in 



the metric of the two-body system ([24| -((26 |) . is different 
from the potential V in Eq. (pOj) . which is the metric 
of a single star. Since these potentials arc written as 
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expansions of powers of 1/c", in the following wc shall 
identify the order of expansion with a superscript index. 
Thus, defines the scalar potential of order in 1/c, 
V^^'^ is the term and so on. 



C. The orbital motion 

Following [42|, we assume that the orbit evolves as a 
slow adiabatic inspiral of a quasi-circular orbit, i.e., the 
energy lost through gravitational waves is balanced by a 
change of the total binding energy E of the system 



dE 
'dt 



(27) 



where E and the GW flux T are expressed in terms of 
the PN variable 



\ 2/3 



(28) 



being ui 



d(j}/dt the orbital frequency. Eq. ([27| yields 
J" 



dx 
IE 



' dE/dx 



(29) 



We neglect the orbital eccentricity because, due to grav- 
itational wave emission, the orbit circularizes well before 
the latest stages of the inspiral which we are studying 
(43| . We use the approach named "Taylor T4 approxi- 
mant" [13, [H] , in which the right-hand side of eq. ([291) is 
expanded to 3.5 PN order including spin terms. We also 
include the effects of the NS tidal deformation on the or- 
bital motion, up to next-to- leading-order Q . The orbital 
phase 0(t) and the orbital frequency ui arc computed by 
numerically integrating the following ODEs 



dx 
dcf) 



dx 



pp 



dx 



tidal 



Gm ■ 



where the point-particle contribution reads 



dx 
'dt 



PP 



5 m 



akX 



fe/2 



(30) 
(31) 

(32) 



k=0 



with the coefficient given in Appendix|Xl and the tidal 
term is given by 



dx 
It 



32miA2 
5m'^ 



12 



[l + ll^ 

L 7Tt 



10 I 

X + 



4421 12263 m2 



28 



28 



+ 



1893 mi 



661 



mi 



1 o 2 



(33) 



where A2 is the Love number of the body 2, and 1 2 
means the same terms but whit the label 1 and 2 ex- 
changed. As wc discuss in Section llll Bl the values of the 
Love number for different stellar models can be computed 



with our approach, and agree with the values obtained 
in the literature 01 • 

The orbital separation ri2 is evaluated through the PN 
expression for 7 = Gm/ri2C^ , which is known up to order 
3 PN, including spin terms [47[, and is found solving the 
equation 



d'y 
Iti 



dx 
'dt 

3.t2 
229 



2x11-- 



5^3/2 ( 



65 
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12 



4x' 

1/2 

^81 



10 8 

1 V I Sf- 

3 9 ' 



2203 41 



2520 192 



(34) 



where S = ^ — ^ia. and the spin variables are defined as 
follows 



G 



c Si 



G'. 



G m 
c 



m2 



771 1 



(35) 
(36) 



Si — {G/c)mfdiSi arc the spin angular momenta of bod- 
ies i = 1,2, with dimensionless spin parameters di and 
unit direction vectors s^. 

It is important to remark that the adiabatic inspiral 
of the orbital motion and the "adiabatic approximation" 
for the Love number, are two different approximations: 
the first assumes that the orbital timescale is much larger 
than that associated to the gravitational wave energy loss 
(orbital adiabatic approximation); the second assumes, as 
mentioned in the Introduction and in Section llll Bl that 
the orbital timescale is much larger than the timescale as- 
sociated to the NS internal dynamics (Love number adi- 
abatic approximation). In this paper, we use the orbital 
adiabatic approximation, but we drop the Love number 
adiabatic approximation. 



D. Post-Newtonian tidal deformations 

Tidal interactions in binary systems have been studied 
by many authors in the framework of general relativity 
(see for instance [H, H^). They are described by the 
equation of geodesic deviation: 



15^ 



, 



(37) 



li'^Vp, and. 



where R°'fj^s Riemann tensor, D / Dt 

in the present case, is the 4— velocity of the star center 
O* and ^" is the separation 4— vector between O* and 
a generic fluid element. By introducing an orthonormal 
tetrad field {e^'j-j} (i = 0, . . . , 3) associated with the frame 
centered in O* , parallel transported along its motion, and 
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such that ej^Q^ = u'', Eq. (|37|) can be cast in the form [5C 



(2) 



-2^(4) 



(38) 



where the ^'^^^ = ejj''^'^, and C'''(]) are the components of 
the relativistic tidal tensor, defined in Eq. ([1]). In the 
affinc approach, Eq. (|38p apphes with replaced by a^. 

In the following subsections, starting from the 3 PN 
metric given in Eqns. ([24)) -p6 l) . we write the explicit ex- 
pression of the parallel transported tetrad, and compute 
the tidal tensor assuming equatorial motion. 



1. The parallel transported tetrad 

The orthonormal tetrad associated to the PN metric 
- (|26|) . satisfies the Fermi- Walker transport equations 
(expressed in terms of the coordinate time t, rather than 
of the proper time of O* ) [lH : 



dt " 



where 



(39) 



(40) 



a'' is the 4— acceleration of O* and = dx^/dt, fi = 0,3 
its coordinate velocity, with v° = c. The tetrad vectors 
are ISlI: 



where 





= 1 + 








c 




c 




= 5i 



V - 



V - 



= — +0(3) 



0(4) 
+ 0(5) 

0(4) 



(41) 



and 



1 



X = ^^xy (42) 
is the angle describing geodesic precession and frame 
dragging, given in terms of the antisymmetric matrix Q 
defined as 

Q(t, to) = / [v X (VF - a) - V X (Fv - 2V)] dt . 

J to 

(43) 

V = {Vi} is the post Newtonian potential, associated 
with the components goi of the PN metric ([25|) . and to is 
an arbitrary integration constant. 

Eqns. (|4T|) reduce to those given in ref. [5l| with the 
identification — jSijV — Xij and Vi 



2. The tidal tensor 



e(t) = 6(4) 
pt pt 



-{x) 
J 



-sinxe|^) -f cosxe^^^) 



Having defined the tetrad field, we have explicitly com- 
puted (with the help of the symbolic manipulation soft- 
ware maple and the package GRTensor) the Riemann and 
the tidal tensors, up to order The general structure 

of the tidal tensor components in terms of the derivatives 
of the PN potentials, is given in Appendix |B] here we 
show, as an example, the component C(^^'^(^^y. 



(x)(;r) 



where l40|, |41| 

^(0) ^ ^ 



2vl 



1^2 



1 o 2 



(45) 



y(3) = 



Gm2 



4riri2 4rirf2 
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2ri 



-2Ge,,kv\Sidk ( - I + 1 2 
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To hereafter we omit the parentheses to indicate the 
tetrad components of the tidal tensor. In Eqns. (|45)) 
1 O 2 means the same term but with the labels 1 and 2 
exchanged; ri = |x — y;^| and rii = (x — yi)/ri, where 
X is the field point and Yi{t) the trajectory of mi (and 
similarly for TO2) ; Vi = dy^{t)/dt is the coordinate veloc- 
ity, ri2 = ri — r2 the relative displacement between the 
two masses, and V12 = vi — V2 the relative velocity, eijk 
is the 3— dimensional antisymmetric Levi-Civita symbol 
with 6123 = Ij find 5*1 is the spin- vector. 

The steps needed to evaluate the tidal tensor at the 
center of the NS (i.e., at location 1 in our conventions) 
are the following: 



1. estimate the various derivatives of the PN poten- 
tials V,Vr, 



2. compute the tidal tensor at the source location 
[Cj]i = Cj(x — >• yi), and apply a regularization 
procedure; 



form 



i^(x,yi,y2) = ^rfA(ni,yi,y2) keZ. (46) 

k 

The regularized value of F at the point 1 is the Hadamard 
part finie, which is the average, with respect to the di- 
rection rii, of the fc = term in the sum (|46|) : 



(^^)i=^^(yi,yi,y2)= / ^^./o(ni,yi,y2) . (47) 

We use this procedure to evaluate {V)i, {Vi)i and their 
derivatives at the source point. We remark that the regu- 
larization procedure should be applied separately on the 
Ricmann tensor and on the orthonormal tetrad. How- 
ever, at the PN order we arc considering, it is perfectly 
equivalent to apply the regularization procedure directly 



to the tidal tensor Cf. 



We now express the point particle positions 2 in the 
system center of mass frame, by the following coordinate 
transformation 1521 



Vi 

2/2 



m2 ^ mi - TO2 „ 

m 771 

mi mi — m2 
1_ jy 7? 

m m 



12 



0(4) 



(48) 



12 



0(4) , 



where 



3. express all quantities in the two-body center of 
mass frame; 



4. switch to the star principal frame defined in Sec- 
tion [llXl 



The second point requires further clarifications. Since the 
PN metric refers to pointlike sources, the PN potentials 
(|^5)) diverge when computed at the source locations x 
Yi and x — > y2. To remove this divergence, we apply the 
Hadamard regularization procedure j40| , which we briefly 
describe. Let -^(x, , y2) be a function depending on the 
field point x and on the two source locations yi,y2, and 
admitting, when x approaches yi, an expansion of the 

I 



V ^ — 

9 



^12 



Gr 



2ri2 



0(4) 



(49) 



Finally, we express Cj in the principal frame using the 
rotation matrix 



T 



cos Tp sin Tp 
— sin tp cos tp 
1 



(50) 



where "0 is defined by Eq. ^ 
dr 



(51) 



The complete form of the tidal tensor c = TCT'^ is given 
by: 



Cxx = - {1 -I- 3cos[2iljt]} + .f^ l [GGml + 5Gmfi + Srl2'miuri2] (1 + 3 cos[2i/';]) 



-I- 6m2ri2 ( ^ +2u^ 0ri2 sin[2V'i] [ 4- J^^'^fj { (j>{rn2 - mi) + (7712 - 5mi)cpcos[2ipi] + ri2 



^m2r-?2(H-cos[2V'i]) + 
(m2 -I- 3mi^ 



ri2 



Cyy = — {1 ^ 3cos[2V'i]} + ^^^4 i [6Gmi + 5Gmfi + 3f-?2n7iz^ri2] (1 - 3cos[2?/>;]) — 6(j}^m2ri2{l - cos[2ipi]) + 
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6m2ri2 + 2u] 0ri2 sin[2i/'i] \ H (I>im2 - rrii) - 4>{m2 - 5mi) cos[2V;;] - ri2 



(m2 + 3mi) 
ri2 



i[2^i] k53) 



„Gm2 5 Gmn Zm\vr^2 3 ;2 

12 ^ ' 12 ^ ' 12 '12 



r 

3Gm2 gjjjj27/);] + ^J,^^ 1 2m2ri2 + 2;/^ 0fi2 cos[2il)i\ - ^Gml + 5Gm/j + 3fi2mi!/ri2 - 2(t?m2ri2^ 



6G(m2 - mi)g| 



(54) 



2r-?2 



X sin[27/'i] )■ < fi2(m2 + 3mi) cos[27/j/] + 0(m2 - 5mi)ri2 sin[2V'zl \ 



(55) 



where the /at; angle ipi = ijj — (f)+x describes the misahgn- 
mcnt between the axis ai and the hne between the two 
objects. In the tidal tensor components the dot indicates 
differentiation with respect to the coordinate time t. In 
the principal frame, the geodesic deviation equation for 
the tidal deformation can be written as 



+ Cijaj = . 



(56) 



It should be stressed that, as noted in 0], if the 
system is in quasi-circular inspiral, the radial motion 
is due only to gravitational back-reaction; consequently 
fi2 ~ (ni2Vi2) ^ and can be neglected. 

3. Comparison with previous expressions of the tidal tensor 

As a first check, we compare the tidal tensor derived in 
[igj for a test particle (i^ — > 0) moving along a geodesic, 
with our tidal tensor. This tensor has been used in the 
literature to study tidal effects in binary systems using 
a quasi-stationary approach [l^ [2^ , or evolving the or- 
bital equations, assuming quasi-circular orbit |3l| . Let 
us consider, as an example, th e c^ x component for a non- 
rotating BH. Eq. (70) of ref. gives 



„sch 



Gm2 



1 



K 



(57) 



where 



C2 c2 i dT 



dt 



and is the radial distance in Schwarzschild coordinates. 
In order to compare Eq. ([57)1 with Eq. ([5^ we need to ex- 
press df^^ in terms of the same radial coordinate adopted 
for the PN expansion 



Gm2 



(59) 



We find (up to l/c' terms) 

G?7l2 



2rh 



(l + 3cos[2V'i]) + ^^{Gm2 



3Gm2 cos [2^/1; 
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i^rJa cos[2ijji 



} -(60) 



This expression coincides with our Eq. (|52p. in the limit 
v ^ and ri2 — 0. 

We would like to make a further remark about the 
difference between the tidal tensor ([57)1 . derived from the 
Schwarzschild metric assuming that mi follows a time- 
like geodesic of the Schwarzschild spacetime, and that 
derived from a two-body post-Newtonian metric. For a 
particle in circular orbit, the constant K given in Eq. 



K 



Gm2 



3Gm2 



(61) 



The former equation diverges for rg — ^ 3Gto2/c^- This 
divergence is present also in the tidal tensor components, 
as shown by Eq. ([57| , and it may affect the evaluation of 
tidal effects even if the distance between the interacting 
bodies is larger than (but close to) rg = 3Gm2/c^. 

Conversely, as stressed in Q , such divergence does not 
appear in the PN equations of motion, and consequently 
the tidal tensor components ([52|l - ([55|l are free of this un- 
physical behaviour. 

On the other hand, as — > our approach loses accu- 
racy, since in the test particle limit the PN expansion is 
poorly convergent 

As a second check we compare our tidal tensor with 
that used in previously derived in up to order 
^ 1 jf? with a completely different approach, based on a 
multipole expansion. Comparing ~G^l (Eq. (2.2) of Q) 
with our tensor Cy, truncated to order ^ l/c^, we find 
that (renaming rri\ <-> m-i) they coincide. 



E. Internal dynamics 

The internal dynamics of the NS is described using 
the Hamiltonian approach in the affine approximation 
[23I . , recently improved to take into account general 
relativistic effects pSj : 



(62) 
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where T-Li describes the NS internal structure, and T-Lt 
describes the tidal interaction. Hj is obtained directly 
from the internal Lagrangian Cj jS]). The tidal Hamilto- 
nian T-Lt is obtained from the tidal Lagrangian (built up 
with the coefficients Cy , Eqns. ([52|) - ([54)) ): 



Ct 



(63) 



where 



\ = M. ■ diag 



Rns 



is the inertia tensor of the star in the principal frame. 

In deriving the dynamical equations from the Hamil- 
tonian ([62ll . we use the PN time coordinate t, which is 
related to the proper time of the star center of mass r by 
the relation dr ~ ^{t)~^dt where the redshift factor 7(t) 
is: 



7(t) = 1 +3 



Gr 



ri2 



,4„2 
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mivri2 



777-2 / 3 2 2 
(4777? + ll?77i 7772 + 1477717(7,2 



5777n 



-f (Ami - 47n/i 



(64) 



It should be mentioned that in previous works, where 
the afhne approach including relativistic corrections was 
used [l^ [2^ [3l| , the contribution of the redshift factor 
was neglected, i.e., it was assumed t ~ r. 

We also remark that Ht — —^t, since Ct does not 
depend on the conjugate momenta. 

The equations of motion for the variables qi = 
{7/), A, ai, 02, 03} and their conjugate momenta pi = 

) P\ ) Pai 1 Pa2 J Pas } are: 



da I 
~dt 
da2 
~dt 

dt 

dpai 
dt 



dpa^ 
dt 



Rns Pai 
lit) M 
Rns Pa2 

Rns Pas 
lit) M 

M 

lit) 



(65) 
(66) 
(67) 



A" 



M al 



ai 



ai 



1 V 3 ~ 
^ M 



(68) 



M 

lit) 



^Ns n 

M a\ 



^yy 



2^Kn + \^R%sA2 

0.2 ^ M. 



02 



(69) 



dpas 

dt 

dX 
dt 
dp\ 
dt 
d4_ 
Itt 

dp^ 

~dr 



M 

lit) 
A 

W) 

1 dC 

jifjdT 

n 

W) 
1 dj 



iv_ 
2m 



-—Rns^s 



E^rg II 
M ol 



= 



M Cxy 



Jit) dT Rns lit) 



where 

A, E 



du 



ia] + u)y^{al+u)ial+u){al + u) 



03 (70) 

(71) 
(72) 
(73) 
(74) 

(75) 



is the NS angular momentum, and C is the conjugate 
momentum associated to A: 



C = [{al + al) A- 2aia2D.] 



J' 



r>2 

^NS 

M 

??2 

^NS 



[[a\ + 03) - 2aia2A] 



(76) 



C can be interpreted as the circulation of the fluid [25|, 
i.e., the line integral of the four- velocity on a closed 
worldline enclosing the system. In absence of viscosity, C 
is a constant of motion. 

It is worth mentioning that, in the spherical configu- 
ration (oi = Rns), the integrals (|75|) can be solved ana- 
lytically, finding Ai = 2/(3i?^g); furthermore, the virial 



theorem ([19)) implies that V = —311. Consequently, in 
the spherical limit the terms in V and 11 in Eqns. (|68p - 
(fzni) cancel: 



2M 



NS 



A,. 



^Ns n 

M af 



= 0. 



(77) 



- sph 



This property is crucial to ensure a stable evolution. In- 
deed, the system (|65|) - (f74l) admits an equilibrium solu- 
tion, for which the star is spherical and non-rotating, 
and the tidal tensor vanishes, only if the property (j77p 
is satisfied. If such solution exists, the tidal deformation 
induced by the interaction is basically a perturbation of 
the equilibrium configuration, and the system of equa- 
tions is well behaved. Conversely, if the cancellation ([77)1 
is not exact, the system becomes unstable, because the 
terms in V and 11 are larger than other terms and the 
equations are non-linear. We remark that the validity of 
Eq. ([77)) is guaranteed in our approach, because the virial 
theorem is satisfied exactly, as discussed in Section Hi AI 



F. Roche lobe and mass shedding 

In the next Section we shall compare the results ob- 
tained by numerically integrating the equations of motion 
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([65)) - (|74p for a a BH-NS coalescence, with those published 
in the literature; we shall evolve the equations up to an 
orbital separation, rshed, at which mass shedding sets in. 

In order to find the value of rshed, we estimate the 
Roche lobe radius of the NS during the inspiral. It defines 
the region surrounding the star where a particle of mass 
Too ^ TOi is bounded to the NS gravitational attrac- 
tion. Following the strategy adopted in [s^ , we estimate 
the three-body potential for masses toq <ti toi < m2 (at 
Newtonian order) for equatorial orbits in the x^y plane: 



U{x,y) 



Gtoi 



Gm.2 

|x-y2l 



1 



Lo^x^ (78) 



where 'yi/2 are the mi/2 displacement vectors, and the 
last term is the centrifugal contribution with 



/ G(toi + m2) 



' 12 



Since U{x,y) takes its maximum, Uri, on the surface 
defining the Roche lobe, we compute numerically U{x^ y), 
finding the Roche lobe on the x — y plane. Mass shed- 
ding starts when the star, which is stretched along the 
direction of the axis ai, touches the Roche lobe. 

We also determine the location of the 3 PN ICO (In- 
nermost Circular Orbit), rjco] to this aim we minimize 
the total binding energy of the BH-NS system, including 
the spin contribution [43]. If r^hed > rjco, the NS is 
disrupted before the merger. 

It should be noted that the affine approach is intrinsi- 
cally non-linear, therefore it takes into account non-linear 
hydrodynamical effects. It also can describe mode oscilla- 
tions (see, for instance, [26j and, in a similar framework, 
54|). However, it does not account for non-linearities in 
the tidal tensor; since these effects are of the order of 
{RNs/f 12)^ 0, which never exceed ^ 10~^, they can be 
safely discarded. 

Higher multipoles (octupole, etc.) of the tidal field are 
also neglected in the affine approach; they are suppressed 
by a factor ^ {Rns /ri2)^ higher PN orders are sup- 

pressed by a factor ^ m/{ri2c'^). Both these quantities 
can become as large as ^ 0.15 as ri2 rshed- Therefore, 
our approach becomes less accurate in the last stages of 
the inspiral. We plan to improve our model, computing 
higher PN orders and higher multipole contributions to 
the tidal field, in future publications. 



III. DYNAMICAL TESTS 

To validate our approach, we have integrated the dy- 
namical equations and (|55 |) - ([711) to simulate 
BH-NS binary coalescences. In order to compare our 
results with the existing literature, we assume that the 
neutron star is irrotational (i.e., we set C = 0), while the 
black hole can rotate. 



A. Fully relativistic simulations 

We compare our results with fully relativistic simu- 
lations from three groups, who have kindly shared the 
required data with us: the Potsdam group [37j (which 
we denote by AEI); the Urbana group (URB); the 
Kyoto/Tokyo group [l3 (KT). All simulations (including 
ours) use the same P = 2 polytropic equation of state. 
The values of the mass ratio q ~ m2/mi ~ Mbh /M^s, 
of the BH spin parameter a ~ a/MsH, and of the NS 
compactness C = M^s/ Rns are: 

1. AEI simulations: g = 5, a = 0, C = 0.1, 0.125, 0.15; 

2. URB simulations: 9 = 3,5 = 0, 0.75, C = 0.145; 

3. KT simulations: g = 2, 3, a = 0, C = 0.145. 

In order to check the validity of our PN formulae, and 
to determine the time offset tof / between our simulations 
and the fully relativistic simulations, preliminary 
check we compare the orbital motion. It is worth remark- 
ing that the time offset is needed to compare our results 
with the fully relativistic simulations. In particular, it is 
needed to compare, for each simulation, the time tshed 
at which our model predicts the onset of mass shedding, 
with the time at which this occurs in the correspond- 
ing fully relativistic simulation. To this aim we follow 
the startegy adopted in [4^: we demand that the PN 
and the numerical gravitational wave frequencies agree 
at some fiducial frequency Wm, defining toff as the time 
for which 4>pNitoff) = w,„. 

We remark that since the frequency of the to = 2 com- 
ponent of the gravitational wave, ilnw = 20, is a gauge 
invariant quantity (see for instance [4y, [ssl . [56j and refer- 
ences therein), it is an appropriate quantity for our com- 
parisons. On the contrary, it is impossible to directly 
compare the radial coordinates, since the gauge used by 
fully relativistic simulations is different from our gauge, 
and furthermore it changes dynamically during the sim- 
ulation jssf . Note that here to is the harmonic index, not 
the total mass as in the rest of the paper. 

The comparison with the URB and KT data is shown 
in Fig. [1] where we plot flaw versus time (both normal- 
ized to the total mass of the binary). Our profiles are 
indicated in Fig. [T] with a dashed line, which ends at the 
onset of mass-shedding. The URB and KT profiles are 
shown by a continuous line. 

As expected [H, [H-il, the PN (and BOB) descrip- 
tion of the inspiral phase is in good agreement with fully 
relativistic simulations. The oscillations in the URB, KT 
curves shown in Fig. [1] are due to the fact that their 
initial data have a residual eccentricity, while our orbits 
are quasi-circular. A comparison with the AEI data gives 
similar results. 

In order to assess the accuracy of our evaluation of 
the onset of mass shedding, we consider the evolution of 
the NS central density, Pc(t), which is a gauge invariant 
quantity. We compare our profiles, with those evaluated 
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by the AEI group for the parameters indicated above. 
The results are shown in Fig. [21 where the AEI profiles 
are plotted with a solid line, and our data with a dashed 
line. In the AEI curves, at some point the central den- 
sity sharply drops down. This signals the transition to 
a new equilibrium configuration and hence the possible 
transfer of mass from the star to the black hole. Deter- 
mining accurately in the numerical simulations when this 
transfer begins, is not trivial (tiny amounts of matter are 
lost from the stellar surface already at large separations) . 
However, it is reasonable to assume that the transfer of 
mass takes place in the transition between the two differ- 
ent values of the central density and therefore in a time 
interval of 0.75 ms for the C=0.1 model or of 0.25 ms for 
the C=0.15 model (clearly, the smaller the compactness 
the slower the transfer process) [s^l- Overall, therefore, 
we can take the decrease of pc to roughly mark the stage 
in which the star fills up the Roche lobe. The dashed line 
for our Pc{t) ends at the onset of mass shedding. This 
occurs just before the steep decrease of the AEI curves, 
for the models with compactness C = 0.1,0.125. For the 
case C ~ 0.15, r/co is reached before mass-shedding sets 
in; therefore, the dashed line in the bottom panel of Fig. [2] 
ends at an earlier time with respect to the sharp drop of 
the solid line. We also note that the values of the cen- 
tral density in our simulations and in AEFs, agree quite 
weU for the models with C = 0.1 and C = 0.125. For 
the model with C = 0.15, the AEI data show an increase 
of the central density at earlier times, probably due to 
spurious numerical effects. 

As a further check, we have compared our estimate of 
the onset of mass shedding with snapshots of the URB 
and KT simulations. In Fig. O we show these snapshots 
at the time t = tshed, which we evaluate for the different 
models. At that time, in the URB and KT simulations 
the NS starts showing a cusp, indicating a mass flow. We 
remark that this kind of comparison should be considered 
as purely qualitative, since the stellar boundary in the 
snapshots corresponds to a threshold value of the stellar 
density, the choice of which is arbitrary. 

We would also like to remark that to compare the re- 
sults of our approach with those of fully relativistic sim- 
ulations is not an easy task. In particular, a "clean" 
comparison of the stellar shape deserves further inves- 
tigations in close interaction with numerical relativity 
groups, and will be considered in future works. 



tions discussed in [s^, HH): 

2 

Qij = ->^Cij = --fc2-R|/sc^j . (79) 

Eq. (I79|) (which is written in the principal frame) is based 
on the "Love number adiabatic approximation" discussed 
in Section III Al C; since it assumes that the timescale 
of the orbital evolution (and then of the tidal tensor 
changes) is much larger then that needed for the star 
to set into a stationary configuration, this assumption 
may not be correct in the latest phase of the inspiral; 
however, it is satisfied when the star and the companion 
are sufficiently far apart. 

In our model, the NS quadrupole moment in the prin- 
cipal frame is 

= ~ a^Sij) (80) 

where = {al+a2+a^)/i, and Cij is given in Eqns. ([52)) - 
(|55p . In order to compare the Love number ^2 predicted 
by our approach with those determined by Hinderer , 
which we denote by fc|^, we have evaluated ^2 using 
Eqns. ([7^ . ([M)) for the same NS models, assuming a 
polytropic equation of state with different values of the 
adiabatic index F and of the compactness C, by setting 
the binary system at the orbital separation of ri2 ~ 180 
km. As shown in Table [J our results agree with those of 
within a few percent. 

We remark that the Love number approach was gener- 
alized in [s^ . [35I , where other Love numbers were intro- 
duced; however, the leading tidal effect is encoded in fc2. 



c 


r 


k2 




0.10 


1.830 


0.0920 


0.0931 


0.15 


1.830 


0.0551 


0.0577 


0.20 


1.830 


0.0297 


0.0327 


0.10 


2.000 


0.1221 


0.1220 


0.15 


2.000 


0.0767 


0.0776 


0.20 


2.000 


0.0444 


0.0459 


0.10 


2.423 


0.1817 


0.1780 


0.15 


2.423 


0.1198 


0.1170 


0.20 


2.423 


0.0737 


0.0721 



TABLE I. The Love number k2, evaluated for different values 
of the NS compactness C, and of the polytropic index F, is 
compared with the values obtained in [3j for the same stellar 
models. 



B. Love number 

A different kind of check can be performed by evalu- 
ating the Love number which, as discussed in the Intro- 
duction, encodes the deformation properties of the star. 
When a weak tidal field induces a deformation on a spher- 
ical star, the star (traceless) quadrupole moment is pro- 
portional to the tidal field [y, |3l (see also the gcneraliza- 



IV. CONCLUDING REMARKS 

In this article we have developed a post-Newtonian- 
affine (PNA) approach which allows to model the tidal 
deformations of a neutron star in compact binary coales- 
cences. To validate the model through a comparison with 
the results of fully relativistic, numerical simulations, we 
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FIG. 1. (Color online) We plot the gravitational wave frequency riow versus time, both normalized to the total mass of the 
binary. The URB, KT data are indicated by a solid line, and our data by a dashed line. The dashed lines stop at ri2 = r^hed, 
where the deformed star touches the Roche lobe. 



have solved the dynamical equations for BH-NS binary 
systems. The tests we have made show a good agreement 
with those results. 

The PNA approach can be useful in many respects. It 
may complement numerical relativity studies of binary 
coalescence because, due to its much lower computational 
cost, it enables to study a large set of models, exploring 
a wide range of parameters. Furthermore, like all semi- 
analytic approaches, it would be helpful to dig the phys- 
ical features of the process out of the numerical artifacts 
which may affect the fully relativistic simulations. 

Since the PNA approach does not assume the "Love 
number adiabatic approximation" , it would allow to test 
the validity domain of this assumption. Indeed, using the 
PNA framework, it would be possible to determine under 
which conditions the NS deformation is characterized by 
a set of constant coefficients, and to find their behaviour, 
if they change during the inspiral. 

Finally, we would like to remind that the production 
of initial data for fully relativistic simulations is a very 
delicate task. Typically, initial data arc plagued by spu- 
rious effects like, for instance, a non-physical eccentricity 
(compact binaries are known to circularize well ahead 
the latest stages of the inspiral); it is difficult to produce 
truly general initial data (for instance, with non-aligned 
spins) . The PNA approach could be used to produce ini- 
tial data for fully relativistic simulations, complementing 
existing initial data solvers. 



These aspects, and the extension of the PNA approach 
to the study of NS-NS coalescences, will be the matter of 
future works. 
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Appendix A: Post-Newtonian expressions for the 
orbital motion 

In this Appendix we write explicitly the post- 
Newtonian coefficients of the Taylor T4 approximant 
eg. pop, for spinning bodies in quasi-circular orbits, with 
spins aligned with the direction of the Newtonian orbital 
angular momentum vector (45j : 
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FIG. 2. (Color online) The NS central density (normalized to the squared total mass of the binary) is plotted, as a function of 
time (normalized to the binary total mass), for the AEI simulations (solid line) and for our simulations (dashed line). 
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Here 7^ is Euler's constant and x = + ^^2, with potentials: 

ai_2 diniensionless spin parameters defined in section lTl CI 



Appendix B: The tidal tensor 



We show the non-vanishing components of the tidal 
tensor Cj up to the Xjc" order, as functions of the PN 
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FIG. 3. (Color online) Snapshots corresponding to t = tnhed which we evaluate integrating our equations. Upper panels: URB 
simulations with C — 0.145, 2 = 0, 0.75, q — 3. Lower panels: KT simulations with C — 0.145, a = 0, q = 2, 3. 
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+ ^ ^{V'^dly - V^dyy ~ dyt)V^^^ + {v^O^y " v'^O^^ " 9rf ) ^ ^ " ^ 9^ y ^ ^ ^ 
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